Documentation
knitr::opts_chunk$set(echo = TRUE)
library(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──## ✓ ggplot2 3.3.2 ✓ purrr 0.3.4
## ✓ tibble 3.0.4 ✓ dplyr 1.0.2
## ✓ tidyr 1.1.2 ✓ stringr 1.4.0
## ✓ readr 1.4.0 ✓ forcats 0.5.0## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(tidycensus)
library(sf)## Linking to GEOS 3.8.1, GDAL 3.1.1, PROJ 6.3.1library(beeswarm)
library(knitr)
library(stargazer)##
## Please cite as:## Hlavac, Marek (2018). stargazer: Well-Formatted Regression and Summary Statistics Tables.## R package version 5.2.2. https://CRAN.R-project.org/package=stargazerlibrary(readxl)
library(mapview)
library(viridis)## Loading required package: viridisLitelibrary(DT)
library(scales)##
## Attaching package: 'scales'## The following object is masked from 'package:viridis':
##
## viridis_pal## The following object is masked from 'package:purrr':
##
## discard## The following object is masked from 'package:readr':
##
## col_factorlibrary(vtable)## Loading required package: kableExtra##
## Attaching package: 'kableExtra'## The following object is masked from 'package:dplyr':
##
## group_rowslibrary(lfe)## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpacklibrary(ggpubr)
mapviewPWS <- function(pwsid,...){
url <- paste0("https://geoconnex.us/ref/pws/",pwsid)
pws <- sf::read_sf(url)
mapview::mapview(pws,...)
}
mapviewOptions(fgb=FALSE)Data Preperation
Find utilities by ownership type
This is the inventory of all utilities that fill out the EAR
inv <- read_excel("Data/EAR/inventory.xlsx",
sheet = "SDWISWaterSystemInventory071620") %>% filter(WaterSystemType=="C")
bounds <- sf::read_sf("https://opendata.arcgis.com/datasets/fbba842bf134497c9d611ad506ec48cc_0.geojson")We have to correct the rates data PWSIDs
rates <- readr::read_csv("Data/OWRS/summary_table_cleanedv2.csv")##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## .default = col_double(),
## utility_name = col_character(),
## pwsid = col_character(),
## effective_date = col_character(),
## bill_frequency = col_character(),
## bill_unit = col_character(),
## bill_type = col_character(),
## tier_starts = col_character(),
## tier_prices = col_character(),
## median_income_category = col_character(),
## bill_frequency_update = col_character()
## )
## ℹ Use `spec()` for the full column specifications.corr <- read_csv("Data/OWRS/correction.csv")##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## pwsid = col_character(),
## utility_name = col_character(),
## `corrected rates PWSID` = col_character()
## )rates <- left_join(rates,corr,by=c("utility_name"))
rates$pwsid.y <- NULL
rates <- rename(rates,pwsid=pwsid.x)
rates$pwsid[which(!is.na(rates$`corrected rates PWSID`))] <- rates$`corrected rates PWSID`[which(!is.na(rates$`corrected rates PWSID`))]Now have to correct rate tier prices
# merge in original bill units
original_units <- read_csv("Data/OWRS/summary_table.csv") %>% select(utility_name,bill_unit) ##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## utility_name = col_character(),
## pwsid = col_character(),
## effective_date = col_date(format = ""),
## bill_frequency = col_character(),
## bill_unit = col_character(),
## usage_ccf = col_double(),
## avg_household_size = col_double(),
## et_amount = col_double(),
## bill_type = col_character(),
## service_charge = col_double(),
## commodity_charge = col_double(),
## bill = col_double(),
## tier_starts = col_character(),
## tier_prices = col_character(),
## median_income_category = col_character(),
## approximate_median_income = col_double()
## )rates <- left_join(rates,original_units,by="utility_name")
#recalculate tier prices that were converted from ccf to kg with appropriate conversion factor dividing by 0.748052kgal/ccf to convert back to $/ccf, and then again to convert to $/kgal. Original calculation was in a wrong unit of $/(ccf/0.748)
rates$Tier_price[which(rates$bill_unit.y=="ccf")]<-rates$Tier_price[which(rates$bill_unit.y=="ccf")]/(0.748052^2)Create Billing data
Now calculate prep to calculate water bills
# set values for non-discretionary water use
gpcd <- 50
days_month <- 365/12
hh_cons_person <- gpcd*days_month
#expand grid of all combinations of utility and household size (1-7)
names <- c(unique(rates$utility_name))
hhsize <- c(1:7)
hhsize <- expand.grid(hhsize,names)
names(hhsize) <- c("hhsize","utility_name")
hhsize$utility_name <- as.character(hhsize$utility_name)
#join to main table
rates <- left_join(rates,hhsize,by="utility_name")
rates$consumption <- (rates$hhsize * hh_cons_person)/1000write function to calculate water bills
rates <- arrange(rates,pwsid,utility_name,hhsize,Tier_number)
rates$Tier_volume[which(!is.na(rates$Tier_price) & is.na(rates$Tier_volume))] <- 10^10
rates <- rates %>%
group_by(utility_name,hhsize) %>%
mutate(next_tier_volume = lead(Tier_volume),
prev_tier_volume = case_when(
Tier_number == 1 ~ 0,
Tier_number > 1 ~ lag(Tier_volume)
),
Tier_width = case_when(
Tier_number == 1 & !is.na(next_tier_volume) ~ Tier_volume,
Tier_number >1 ~ Tier_volume - prev_tier_volume,
is.na(next_tier_volume) ~ 10^10,
),
vol_in_tier = case_when(
bill_type == "Uniform" & Tier_number ==1 ~ consumption,
consumption <= prev_tier_volume ~ 0,
consumption > Tier_volume ~ Tier_width,
consumption > prev_tier_volume & consumption <= Tier_volume ~ consumption - prev_tier_volume#,
# consumption > Tier_volume ~ Tier_width,
#consumption >= Tier_volume & is.na(next_tier_volume) ~ consumption - Tier_volume#,
#consumption >= Tier_volume & is.na(next_tier_volume) ~ consumption - Tier_volume
)
)
rates$vol_revenue <- NA
rates$vol_revenue = rates$Tier_price * rates$vol_in_tier
rates$vol_revenue[which(rates$bill_type=="1.49*usage_ccf")] <- 1.49*rates$consumption[which(rates$bill_type=="1.49*usage_ccf")]*0.748052
rates$vol_revenue[which(rates$bill_type=="3.04*usage_ccf*1.333")] <- 3.04*rates$consumption[which(rates$bill_type=="3.04*usage_ccf*1.333")]*1.333*0.748052
bills <- rates %>% group_by(pwsid,
utility_name,
bill_type,
hhsize,
consumption,
service_charge) %>%
summarise(volumetric_charge = sum(vol_revenue,na.rm=TRUE))## `summarise()` regrouping output by 'pwsid', 'utility_name', 'bill_type', 'hhsize', 'consumption' (override with `.groups` argument)bills$bill <- bills$service_charge + bills$volumetric_chargeNow we need to categorize them all by ownership type. To Recover ownership type, we read in the latest EAR
# ear13 <- read_tsv("Data/EAR/earsurveyresults_2013ry.zip",col_types="ccdcccdddc")
# ear14 <- read_tsv("Data/EAR/earsurveyresults_2014ry.zip",col_types="ccdcccdddc")
# ear15 <- read_tsv("Data/EAR/earsurveyresults_2015ry.zip",col_types="ccdcccdddc")
# ear16 <- read_tsv("Data/EAR/earsurveyresults_2016ry.zip",col_types="ccdcccdddc")
# ear17 <- read_tsv("Data/EAR/earsurveyresults_2017ry.zip",col_types="ccdcccdddc")
# ear18 <- read_tsv("Data/EAR/earsurveyresults_2018ry.zip",col_types="ccdcccdddc")
# ear19 <- read_tsv("Data/EAR/earsurveyresults_2019ry.zip",col_names=c("PWSID",
# "Survey",
# "Year",
# "SectionName",
# "QuestionName",
# "QuestionResults",
# "SurveyID",
# "SectionID",
# "QuestionID",
# "Pk"), col_types="ccdcccdddc")
# ear <- bind_rows(ear13,ear14,ear15,ear16,ear17,ear18,ear19)
# save(ear,file="Data/ear.rds")
load("Data/ear.rds")
ownership <- ear %>% filter(QuestionName %in% c("Water System Ownership",
"IsWholesaler")) %>%
select(PWSID,
Year,
QuestionName,
QuestionResults,
QuestionID,
Pk) %>%
pivot_wider(id_cols = c(PWSID, Year),
names_from = QuestionName,
values_from = QuestionResults,
names_repair="unique") %>% unnest(cols=c("Water System Ownership","IsWholesaler") ) %>%
group_by(PWSID) %>%
mutate(maxYear = max(Year)) %>%
ungroup() %>%
filter(`Water System Ownership` %in% c(
"Privately owned business (non-community)",
"Privately owned Mutual Water Company or Association",
"Privately owned, non-PUC-regulated (Community Water System)",
"Local Government",
"Privately owned, PUC-regulated, for profit water company"
)) ## Warning: Values are not uniquely identified; output will contain list-cols.
## * Use `values_fn = list` to suppress this warning.
## * Use `values_fn = length` to identify where the duplicates arise
## * Use `values_fn = {summary_fun}` to summarise duplicatesCensus data
Get all the variables
options(tigris_use_cache = TRUE)
vars <- load_variables(2019, "acs5", cache = TRUE)
census_vars2 <- c(median_income="B19013_001",
pop_in_hh_with_income_assistance="B09010_002",
pop_in_hh_with_no_income_assistance="B09010_008",
hh_inc_5k="B19001_002",
hh_inc_12.5k="B19001_003",
hh_inc_17.5k="B19001_004",
hh_inc_22.5k="B19001_005",
hh_inc_27.5k="B19001_006",
hh_inc_32.5k="B19001_007",
hh_inc_37.5k="B19001_008",
hh_inc_42.5k="B19001_009",
hh_inc_47.5k="B19001_010",
hh_inc_55k="B19001_011",
hh_inc_67.5k="B19001_012",
hh_inc_87.5k="B19001_013",
hh_inc_112.5k="B19001_014",
hh_inc_137.5k="B19001_015",
hh_inc_175k="B19001_016",
hh_inc_gr200k="B19001_017",
hh_income_count = "B19001_001")Lets get all the Municipal boundaries
st <- get_acs(year=2019,geography = "state",
variables = "B01003_001", geometry = TRUE) %>% filter(NAME=="California")## Getting data from the 2015-2019 5-year ACSpl <- get_acs(year=2019,geography = "place",
variables = "B01003_001", geometry = TRUE, state= "CA") ## Getting data from the 2015-2019 5-year ACSvars <- load_variables(2019, "acs5", cache = TRUE)Let’s get all the Census Block Groups
# bg <- get_acs(year=2019,geography = "block group",
# variables = census_vars2, geometry = TRUE, state= "CA")
# save(bg, file="Data/census_blockgroups.rds")
load("Data/census_blockgroups.rds")Let’s get all the Urbanized Areas and combine them into super areas.
ua <- get_acs(year=2019,geography = "urban area",
variables = c("B01003_001","B19013_001"), geometry = TRUE, output="wide")## Getting data from the 2015-2019 5-year ACS## Warning: `funs()` is deprecated as of dplyr 0.8.0.
## Please use a list of either functions or lambdas:
##
## # Simple named list:
## list(mean = mean, median = median)
##
## # Auto named with `tibble::lst()`:
## tibble::lst(mean, median)
##
## # Using lambdas
## list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.ua<-rename(ua,Pop=B01003_001E)
ua <- ua[st,]## although coordinates are longitude/latitude, st_intersects assumes that they are planarDescriptive stats
Now we categorize all utilities in the SDWIS Inventory by type
own <- ownership %>% group_by(PWSID,`Water System Ownership`) %>% add_tally() %>% ungroup() %>% group_by(PWSID) %>% add_tally() %>% ungroup()## Storing counts in `nn`, as `n` already present in input
## ℹ Use `name = "new_name"` to pick a new name.own1 <- filter(own,n == nn)
own1 <- distinct(own,PWSID,.keep_all=TRUE)
own2 <- filter(own, n != nn)
own2 <- own2 %>%
distinct(PWSID, `Water System Ownership`, .keep_all=TRUE) %>%
group_by(PWSID) %>%
mutate(maxYear = max(Year)) %>%
ungroup() %>%
filter(Year==maxYear) %>%
group_by(PWSID) %>%
add_tally()## Storing counts in `nnn`, as `n` already present in input
## ℹ Use `name = "new_name"` to pick a new name.own1 <- select(own1, PWSID, `Water System Ownership`,Year)
own2 <- select(own2, PWSID, `Water System Ownership`,Year)
own <- bind_rows(own1,own2)
own <- own %>% group_by(PWSID) %>% mutate(maxYear=max(Year)) %>% ungroup() %>% filter(Year == maxYear)
rm(own1,own2)
d <- left_join(inv,own, by="PWSID")
d <- select(d, PWSID, PWS_Name, County, Pop, ServConn, `Water System Ownership`)
d$log10PopCat <- floor(log10(d$Pop))some distributive plots
d <- filter(d,`Water System Ownership` != "Privately owned business (non-community)")
x <- data.frame(table(d$log10PopCat,d$`Water System Ownership`))
x <- x%>%pivot_wider(id_cols = Var1,names_from = Var2, values_from=Freq)
datatable(x,extensions="Buttons",options = list(
dom = 'Bfrtip',
buttons = c('copy', 'csv')))y <- d %>%
mutate(TotPop = sum(Pop,na.rm=TRUE),
PropPop = Pop/TotPop) %>%
group_by(log10PopCat,`Water System Ownership`) %>%
summarise(Population = sum(Pop,na.rm=TRUE),
Pop_percent = 100*sum(PropPop,na.rm=TRUE)) %>%
pivot_wider(id_cols = `log10PopCat`,
names_from = `Water System Ownership`,
values_from = c(Population,Pop_percent))## `summarise()` regrouping output by 'log10PopCat' (override with `.groups` argument)datatable(y,extensions="Buttons",options = list(
dom = 'Bfrtip',
buttons = c('copy', 'csv')))p<- ggplot(d, aes(x=log10(Pop), color=`Water System Ownership`), size=4) +
geom_density(lwd=2) +
scale_color_viridis_d() +
theme_dark() +
theme(legend.position=c(0.69, 0.85)) +
xlab("log(base 10) of Service Population")
p## Warning: Removed 36 rows containing non-finite values (stat_density).
Spatial demonstration figure. Show utility boundaries, municipal boundaries, and censusblock groups
b <- bounds %>% select(SABL_PWSID) %>% rename(PWSID=SABL_PWSID)
d.b <- d %>% inner_join(b,by="PWSID") %>% st_as_sf()
pal <- magma(n = length(unique(pl$NAME)), direction = -1)
bkfield <- st_transform(filter(ua,GEOID=="04681"),4326)
d1 <- d.b[bkfield,,op=st_intersects]## although coordinates are longitude/latitude, st_intersects assumes that they are planarpl1 <- st_transform(pl,4326)[bkfield,,op=st_intersects]## although coordinates are longitude/latitude, st_intersects assumes that they are planarm <- mapview(pl1,
stroke=TRUE,
color="red",
alpha.regions=0,
layer.name = "Municipal Boundaries",
lwd=2,legend=TRUE,col.regions="red") +
mapview(d1,
zcol="Water System Ownership",
alpha=0.8, layer.name = "Water System Boundaries")
#mapshot(m,"Figures/Map1.html")
md.x <- st_buffer(st_transform(d.b,3310),dist=0)Bakersfield example
bkfield <- filter(ua,GEOID=="04681")
bg.bounds <- st_transform(bg.bounds,4326)
bkfield <- st_transform(bkfield,4326)
bkfield <- bg.bounds[bkfield,]## although coordinates are longitude/latitude, st_intersects assumes that they are planarbk.bound <- d.b[bkfield,,op=st_intersects]## although coordinates are longitude/latitude, st_intersects assumes that they are planarm2<-mapview(bkfield,zcol="median_income",layer.name="Census Block Group - Median HH Income") + mapview(bk.bound, alpha=1, layer.name = "Water System Boundaries", stroke=TRUE,color="black",col.regions="black", alpha.regions=0, lwd=4,legend=TRUE)
m2#mapshot(m2,"Figures/Map2.html")Calculate affordability
bg.bounds <- st_transform(bg.bounds,3310)
bg.bounds$Area_BlockGroupPWS <- st_area(bg.bounds)
bg.bounds$Percent_BlockGroupinPWS <- bg.bounds$Area_BlockGroupPWS/bg.bounds$AreaBlockGroup
pws.income <- bg.bounds %>% st_drop_geometry() %>%
left_join(bg,by="GEOID") %>%
filter(variable %in% c('hh_inc_5k',
'hh_inc_12.5k',
'hh_inc_17.5k',
'hh_inc_22.5k',
'hh_inc_27.5k',
'hh_inc_32.5k',
'hh_inc_37.5k',
'hh_inc_42.5k',
'hh_inc_47.5k',
'hh_inc_55k',
'hh_inc_67.5k',
'hh_inc_87.5k',
'hh_inc_112.5k',
'hh_inc_137.5k',
'hh_inc_175k',
'hh_inc_gr200k')) %>%
group_by(PWSID,variable) %>%
summarise(estimate=sum(estimate*Percent_BlockGroupinPWS)) %>%
ungroup() %>%
mutate(var2=case_when(variable=="hh_inc_5k" ~ 1,
variable=="hh_inc_12.5k" ~ 2,
variable=="hh_inc_17.5k" ~ 3,
variable=="hh_inc_22.5k" ~ 4,
variable=="hh_inc_27.5k" ~ 5,
variable=="hh_inc_32.5k" ~ 6,
variable=="hh_inc_37.5k" ~ 7,
variable=="hh_inc_42.5k" ~ 8,
variable=="hh_inc_47.5k" ~ 9,
variable=="hh_inc_55k" ~ 10,
variable=="hh_inc_67.5k" ~ 11,
variable=="hh_inc_87.5k" ~ 12,
variable=="hh_inc_112.5k" ~ 13,
variable=="hh_inc_137.5k" ~ 14,
variable=="hh_inc_175k" ~ 15,
variable=="hh_inc_gr200k" ~ 16,)) %>%
arrange(PWSID,var2) %>%
group_by(PWSID) %>%
mutate(cum_inc_sum = cumsum(estimate),
percentile_inc = 100*cum_inc_sum/sum(estimate)) %>%
ungroup() %>%
mutate(perc_dist20=abs(as.numeric(percentile_inc)-20),
perc_dist50=abs(as.numeric(percentile_inc)-50)) %>%
group_by(PWSID) %>%
filter(perc_dist20==min(perc_dist20) |
perc_dist50==min(perc_dist50)) %>%
mutate(percentile=case_when(var2==min(var2)~"20th Percentile",
var2==max(var2)~"50th Percentile"),
income_perc=case_when(var2==1 ~ 5000,
var2==2 ~ 12500,
var2==3 ~ 17500,
var2==4 ~ 22500,
var2==5 ~ 27500,
var2==6 ~ 32500,
var2==7 ~ 37500,
var2==8 ~ 42500,
var2==9 ~ 47500,
var2==10 ~ 55000,
var2==11 ~ 67500,
var2==12 ~ 87500,
var2==13 ~ 112500,
var2==14 ~ 137500,
var2==15 ~ 175000,
var2==16 ~ 225000)
) %>%
pivot_wider(id_cols=PWSID,names_from=percentile,values_from=income_perc)## `summarise()` regrouping output by 'PWSID' (override with `.groups` argument)## Warning: Values are not uniquely identified; output will contain list-cols.
## * Use `values_fn = list` to suppress this warning.
## * Use `values_fn = length` to identify where the duplicates arise
## * Use `values_fn = {summary_fun}` to summarise duplicatespws.income2 <-unnest(pws.income, cols = c(`20th Percentile`, `50th Percentile`, `NA`)) %>%
select(PWSID,`20th Percentile`,`50th Percentile`) %>% distinct(.keep_all=TRUE)
aff <- left_join(bills,pws.income2,by=c("pwsid"="PWSID"))
aff$AR20 <- 12*100*aff$bill/aff$`20th Percentile`
aff$MHI <- 12*100*aff$bill/aff$`50th Percentile`
aff <- rename(aff,PWSID=pwsid)
aff <- left_join(aff,d.x,by="PWSID")
aff <- st_as_sf(aff)
aff <- distinct(aff,PWSID,`Water System Ownership`,hhsize,.keep_all=TRUE)Urban-Rural-Income
d.ua <- st_join(d.x,st_transform(ua,3310), join=st_intersects, largest=TRUE)## Warning: attribute variables are assumed to be spatially constant throughout all
## geometriesd.ua$Pop.y[which(is.na(d.ua$Pop.y))]<-10
d.ua <- distinct(d.ua,PWSID,.keep_all=TRUE)
p<-ggplot(d.ua,aes(x=Pop.y,y=Pop.x, color=`Water System Ownership`)) +
geom_point() +
scale_x_log10() + scale_y_log10() + scale_color_viridis_d() + theme_dark() +
xlab("Population of Urbanized Area") + ylab("Water System Service Population") +
theme(legend.position=c(0.4,0.9))
png("Figures/Fig2.png")
p## Warning: Transformation introduced infinite values in continuous y-axisdev.off()## quartz_off_screen
## 2p## Warning: Transformation introduced infinite values in continuous y-axis
inc <- left_join(pws.income2,d.ua,by="PWSID")
# p<-ggplot(inc,aes(x=Pop.y,y=`50th Percentile`, color=`Water System Ownership`, size=Pop.x/1000)) + scale_size_binned("Service Pop. (1000s)", trans="identity", breaks=c(1,10,100,1000,10000)) +
# geom_point() +
# scale_x_log10() + scale_y_log10() + scale_color_viridis_d(labels = function(x) str_wrap(x, width=13)) + theme_dark() +
# xlab("Population of Urbanized Area") + ylab("Median Household Income within Service Area")
inc$cat <- "1. Pop. <1,000 (n=1753)"
inc$cat[which(inc$Pop.x>=1000)] <- "2. Pop. 1,000 - 10,000 (n=456)"
inc$cat[which(inc$Pop.x>=10000)] <- "3. Pop. 10,000 - 100,000 (n=326)"
inc$cat[which(inc$Pop.x>=100000)] <- "4. Pop. >= 100,000 (n=82)"
p<- ggplot(inc, aes(x=`50th Percentile`, color=`Water System Ownership`), size=4) +
geom_density(lwd=1.24, adjust=2) +
scale_color_viridis_d(labels = function(x) str_wrap(x, width=13)) +
theme_dark() + scale_x_log10() +
xlab("Median HH Income") + facet_wrap(vars(cat))
png("Figures/Fig3.png")
p## Warning: Removed 6 rows containing non-finite values (stat_density).## Warning: Groups with fewer than two data points have been dropped.## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Infdev.off()## quartz_off_screen
## 2p## Warning: Removed 6 rows containing non-finite values (stat_density).## Warning: Groups with fewer than two data points have been dropped.## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning -
## Inf
p<- ggplot(aff, aes(x=`50th Percentile`, color=`Water System Ownership`), size=4) +
geom_density(lwd=2, adjust=2) +
scale_color_viridis_d(labels = function(x) str_wrap(x, width=13)) +
theme_dark() + scale_x_log10() +
xlab("Median HH Income")
png("Figures/Fig3b.png")
p## Warning: Removed 147 rows containing non-finite values (stat_density).dev.off()## quartz_off_screen
## 2p## Warning: Removed 147 rows containing non-finite values (stat_density).
Big map (urbanized areas and utilities)
m3 <- mapview(d.x,zcol="Water System Ownership", alpha=0.8, layer.name = "All Water System Boundaries") +
mapview(filter(aff,hhsize==4),stroke=TRUE,color="red",alpha.regions=0, layer.name = "Systems with price information", alpha = 0.5, col.regions="red")
#mapshot(m3,"Figures/Map3.html")
m3